Concentration fluctuations and phase transitions in coupled modulated bilayers 
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We consider the formation of finite-size domains in lipid bilayers consisting of saturated and hybrid 
lipids. First, we describe a monolayer model that includes a coupling between a compositional 
scalar field and a two-dimensional vectorial order-parameter. Such a coupling yields an effective 
two-dimensional microemulsion free-energy for the lipid monolayer, and its characteristic length 
of compositional modulations can be considered as the origin of finite-size domains in biological 
membranes. Next, we consider a coupled bilayer composed of two modulated monolayers, and discuss 
the static and dynamic properties of concentration fiuctuations above the transition temperature. 
We also investigate the micro-phase separation below the transition temperature, and compare the 
micro-phase separated structures with statics and dynamics of concentration fiuctuations above the 
transition. 
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I. INTRODUCTION 

Biomembranes are two-dimensional (2D) fluids that 
separate the inner and outer environment of organelles in 
biological cells. Naturally occurring biomembranes con- 
sist typically of numerous lipid species, sterols, sugars 
and membrane proteins. According to the "lipid raft" 
hypothesis 1], some of the lipid components and/or pro- 
teins are incorporated into finite-size domains, which play 
an important role on cellular functions such as signal 
transduction processes. Recent experiments suggest that 
lipid rafts are nothing but dynamical molecular assem- 
blies of 20 nm in size with finite lifetime in the order of 
10-20 ms 0]. 

Being motivated by the raft hypothesis, a large num- 
ber of investigations |3l-[l4| have been conducted to reveal 
the properties of artificial membranes consisting of lipid 
mixtures and cholesterol. Below the miscibility transition 
temperature, formation of micron-size domains were ob- 
served using fluorescent microscopy [^-IS]- In some cases, 
rather than a macroscopic phase separation, domains 
with distinct size in the micrometer range have been re- 
ported 0^1. For example, various types of modulated 
(stripe or hexagonal) patterns have been found for multi- 
component lipid and cholesterol mixtures [loj . Above the 
miscibility transition temperature, even multicomponent 
membranes do not phase separate, and their concentra- 
tion fluctuations around the homogeneous state can be 
investigated [ll| - [l3| . in particular, close to the critical 
point, Tc- Furthermore, it is interesting to note that 
critical concentration fluctuations have been observed in 
membranes extracted from living cells 

One of the main reasons that initiated the notion of 
lipid rafts in biomembranes is the existence of finite-size 
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domains, rather than domains resulting from a macro- 
scopic phase separation jj^]. Assuming that membranes 
are in equilibrium, the same question can be phrased in 
terms of whether the biomembrane state is above or be- 
low the miscibility transition temperature. In the high 
temperature one-phase region, the only relevant length 
scale is the correlation length associated with concentra- 
tion fluctuations, and it diverges at the critical point, 
Tc [I6j. Below this temperature, there should be a phys- 
ical mechanism suppressing domain coarsening in order 
to explain the existence of finite-size domains in equilib- 
rium. 

Yet another characteristic feature of biomembranes is 
that the lipid composition of the two leaflets (monolay- 
ers) constituting the bilayer is not the same 17|. More- 
over, such asymmetric monolayers are not independent 
but are coupled to one another. This was confirmed ex- 
perimentally by investigating the phase separation of bi- 
layers with different monolayer lipid compositions J_8J , or 
seen in simulations [l^ . One mechanism that leads to the 
coupling between the two leaflets is the lipid chain inter- 
digitation occurring at the mid-plane of the bilayer (20l - 
|22|, which may affect the domain size in asymmetric bi- 
layers. 

One possibility to account for such finite-size domains 
in lipid mixtures is to consider the special role played by 
"hybrid lipids" [23l - [26| . These lipids have one saturated 
hydrocarbon chain and another unsaturated one. They 
tend to be localized at boundaries of 2D domains and act 
as a line-active agent. Hence, a ternary lipid mixture con- 
sisting of saturated lipids with two saturated chains, un- 
saturated lipids with two unsaturated chains, and hybrid 
lipids can be regarded as a "2D microemulsion" , because 
it is analogous to the microemulsion phase consisting of 
oil, water and surfactant [l^ [13] . For microemulsions it 
is known that there is another length scale, in addition 
to the correlation length mentioned above, related to the 
size of water and oil micro-domains 12 8i . 
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Using the microemulsion analogy, we suggest that an 
additional length scale can be responsible for the finite- 
size structures in biomembranes. We propose a model 
for bilayers consisting of two coupled monolayers that are 
both in a 2D microemulsion state, manifesting modulated 
phases. First, we describe a monolayer model that in- 
cludes a coupling between the lipid composition and a 2D 
vectorial order-parameter. The model accounts for the 
line active nature of hybrid lipids. Next, we consider a 
bilayer in which two such monolayers are coupled through 
an inter-leaflet interaction. For the coupled bilayer, we 
discuss the static and dynamic properties of concentra- 
tion fluctuations above Tc, and investigate their micro- 
phase separation below Tc [2§|. Intermediate structures 
arise when two competing structures have different char- 
acteristic length scales. One of our important conclusions 
is that the micro-phase separated structures below Tc re- 
flect the static and dynamic properties of concentration 
fluctuations above Tc. 

This paper is organized as follows. In Sec. HIl after 
describing how the microemulsion state is obtained for 
monolayers, a model for coupled modulated bilayers is 
presented. In Sec. IIIIl we show the results of static and 
dynamic structure factors for the coupled bilayers above 
Tc. In Sec. IIVI we describe some results for competing 
micro-phase separation in coupled bilayers below Tc. Fi- 
nally, in Sec. |Vl further discussion and flnal remarks are 
presented. 



II. MODEL 
A. Modulated lipid monolayers 

We consider a 2D microemulsion formation in a mono- 
layer consisting of two types of lipids: a saturated lipid 
(denoted by "S") and a hybrid one ("H"). As shown 
schematically in Fig.[lja), the saturated lipid has two sat- 
urated chains, whereas the hybrid lipid has one saturated 
chain and an unsaturated one. It is well known that 
the resulting liquid-ordered (Lq) and liquid-disordered 
(Ld) phases are rich in saturated and hybrid lipids, re- 
spectively, as is depicted in Fig. Hfc). In the experi- 
ments mentioned previously, cholesterol is usually added 
as a third component, and is known to affect the area 
per headgroup of lipids [s^, HH. However, since choles- 
terol has a strong preference for the Lq phase and affects 
mainly the saturated lipid, we do not consider cholesterol 
explicitly in our model, and neglect its presence hereafter. 

The local area fraction of saturated and hybrid lipids 
are defined as (psi^) and 0H(r), respectively, where r = 
(x, y) is a 2D vector. Under the incompressibility condi- 
tion, i)f>s + '/'H = 1, and the only relevant order-parameter 
is the difference between the two lipid compositions, 
4' — 4's~ ■ The phenomenological free-energy that de- 
scribes the 2D phase-separation between S and H lipids 
is given by a Ginzburg-Landau expansion in terms of the 
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FIG. 1. (a) Saturated lipids (S) with two saturated hydro- 
carbon chains, and hybrid lipids (H) with one saturated and 
another unsaturated chain. For each hybrid lipid, we define 
a 2D vector b pointing from its unsaturated tail toward the 
saturated one. (b) The vector m denotes a coarse-grained 2D 
vectorial order-parameter over the microscopic b vectors, (c) 
A binary lipid monolayer consisting of saturated and hybrid 
lipids. At the interface between the Lo and Ld phases, hybrid 
lipids orient their saturated chains toward the Lo phase. The 
distribution of the m(r) vector field is shown by arrows of 
variable length. 



order-parameter 
Fs[4>] = 



+ - /i0 



(1) 



Here cr > is related to the line tension between mono- 
layer domains, r oc T — Tc is the reduced temperature 
with respect to the critical point Tc, and ^ is the chem- 
ical potential that regulates the average <p value in the 
monolayer. Without loss of generality, the coefficient of 
the quartic term is set to be a positive constant by an 
appropriate rescaling of the position variable, r. 

Next we discuss the role of the hybrid lipid and its 
effect on the phase separation. We define a 2D lat- 
eral vector b pointing from the unsaturated tail of the 
hybrid lipid toward its saturated one, as depicted in 
Fig. [TJa). We then introduce a coarse-grained 2D vec- 
torial order-parameter field m(r) = (ma;(r), my(r)), as 
shown in Fig. [ijb) . This vector is the spatial average of 
the b vectors over areas large as compared with molecu- 
lar size, but still small enough as compared with macro- 
scopic scales. The phenomenological free-energy associ- 
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ated with the vectorial field m (and due to the hybrid 
lipids) can be written as an expansion up to quadratic 
order in m 13211: 



Fh[m] = 



d^r 



f(V-m)2 



— m 

2 



(2) 



The coefficient K is the 2D elastic constant, while a is 
taken to be positive so that m = is the stable homoge- 
neous state. An additional term (V x m)^ is allowed by 
symmetry, but is not considered at present for simplicity 
sake. 

At the Lo/Ld interface, hybrid lipids orient their sat- 
urated and unsaturated chains toward the Lq and Ld 
phases, respectively, thereby reducing the chain mis- 
match [2^ [2^ as shown in Fig. [TJc). In fact, the line 
activity of hybrid lipid was predicted [1^ [2^ to be more 
pronounced for binary mixtures of saturated and hybrid 
lipids as compared to ternary mixtures. Within our phe- 
nomenological approach, the role of hybrid lipids can be 
represented by a coupling between the lateral variation 
of 4) and the vectorial field m. To lowest orders, the 
coupling term in the free energy is given by 



-T j d^rm- (V0), 



(3) 



where F is a positive coupling constant because the vec- 
tor m tends to orient towards the Lq domains. The 
total monolayer free-energy is given here by the sum 
of the three terms introduced in Eqs. ([J)-©: — 
Fs -\- Fh -|- Fc. Notice that our model is valid in the weak 
segregation limit (close to T^), because slow spatial vari- 
ation of the order parameters is intrinsically assumed. 
We further note that a similar free energy functional us- 
ing a vectorial order-parameter was proposed for 3D mi- 
croemulsions "s^. 

The total free energy can be conveniently expressed 
using the 2D Fourier transform of 



(q) = j dV0(r)e-*'i- 



(4) 



where q — {qx,qy) is a 2D wavevector, and similarly 
m(q) is the Fourier transform of m(r). Minimizing 
with respect to m, we obtain its optimum value as 



m(q) 



iFq 



a + Kq^ 



(5) 



where q = |q|. By substituting back Eq. ^ into the 
monolayer free energy Fm, the minimized free energy is 
expressed as 



F^[4] 



d^q 1 



(27r)2 2 



2„2 



T + aq 



a + Kq^ 



(q)'^(-q) 



1 



-0 



(6) 



where the Fourier transform has been used only for the 
second-order terms in 0. Expanding the effective binary 
interaction for small we obtain 



Fr^m 



A?q 
+ / dV 



2Bq4 

1 



- 2Aq^ 



(q)'/'(-q) 



where two new parameters are defined 



B 



1^ 



A = 



1 fv 



(7) 



(8) 



When the coupling constant F is small enough, i.e., 
F^/a < o- (or A < 0), the minimum of the Gaussian term 
(first integral) in Eq. ([7]) occurs at g* =0, and is a signa- 
ture of a macroscopic phase separation (for r < 0). On 
the other hand, when F is large enough, i.e., F^/a > a (or 
A > 0), the minimum occurs at a non-zero wavenumber 
q* 7^ given by 




a(l - acr/F2) 



2K 



(9) 



indicating a potential micro-phase separation with q* > 
modulation. 

Taking the inverse Fourier transform, Eq. ([7]) can be 
expressed in position space as 



F„,[</)] = 1 d^ 



+ _ 



(10) 



This is the so-called "2D microemulsion" free-energy for 
a monolayer composed of a binary lipid mixture. When 
A > 0, the negative gradient-squared term favors spa- 
tial modulations, while the positive Laplacian squared 
term with B > Q suppresses modulations. As mentioned 
above, Yamamoto et al. [2^ [11] showed that the effec- 
tive line tension between domains becomes negative for 
membranes consisting of saturated and hybrid lipids. 

More generally, models based on equations similar to 
Eq. (fTO|) have been used successfully in the past to de- 
scribe modulated phases [34| arising in a variety of differ- 
ent biophysical and chemical systems such as Langmuir 
films 35 1, lipid membranes [H, [s^l and diblock copoly- 
mers 



B. Coupled lipid bilayers 

We consider two coupled modulated monolayers form- 
ing a bilayer as shown in Fig. [2j Each monolayer is a 
binary mixture of saturated and hybrid lipids. We de- 
fine two local order-parameters for the two monolayers: 
</* = 0s ~ "/"H and tp = — ipYi, depicted in Fig. [21 
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The coarse-grained free-energy functional is then written 
as [291: 



III. CONCENTRATION FLUCTUATIONS 
ABOVE Tc 



(11) 



The first five terms depend only on cj) and its derivatives 
and describe the upper monolayer in Fig.[3]and its possi- 
ble modulations. These are the same terms as in Eq. PH)) 
for the single monolayer case. Similarly, the latter five V' 
terms describe the lower monolayer, where and are 
the two reduced temperatures, while /i^ and /i^ are the 
corresponding chemical potentials. The last term, Acjjip, 
represents the coupling between the two leaflets with a 
coupling constant A, 



FIG. 2. Schematic illustration of two coupled modulated 
monolayers forming a bilayer membrane. Each monolayer is 
composed of a binary S/H lipid mixture, which can have a 
lateral composition modulation. The relative composition of 
S and H lipids in the upper and lower leaflets are defined by (p 
and ip, respectively. The lipid tails interact across the bilayer 
mid-plane. 



We would like to comment on the physical origin of 
the proposed coupling term, — At/)-!/), in Eq. ([Tl]) . Note 
that this quadratic term is invariant under exchange of 
^ O '0. When A > 0, this term can be obtained from 
a term such as (0 — V')^, which represents a local en- 
ergy penalty when the upper and lower monolayers have 
different compositions [20, 21]. For mixed lipid bilayers, 
such a coupling may result from the conformational con- 
finement of the lipid chains, and hence would have an 
entropic origin [20|. By estimating the degree of lipid 
chain interdigitation, the order of magnitude of the cou- 
pling parameter A can be evaluated [22j. In general, A 
can also be neg ative depending on the specific coupling 
mechanism [3^. However, as the sign of A does not af- 
fect our result, it is sufficient to consider the A > case. 
The other possible higher-order coupling terms which are 
allowed by symmetry are (p'^ip'^, ^^'^ (fiip. However 
these terms do not affect the properties of concentration 
fluctuations in any essential way, and will not be consid- 
ered hereafter. 



Using the bilayer free-energy, Eq. ([TT]), we obtain the 
static and dynamic structure factors, which describe the 
properties of concentration fluctuations for two coupled 
monolayers in their respective one-phase region (above 
Tc). We shall closely follow the formulation of Ref. [iO] in 
which the coupled macro-phase separation was discussed 
for bilayers. 



A. Static structure factor 

The spatially varying 0(r) and V'(r) can be written as 
0(r) = -I- 54i{y) and V'(r) = V'o + 5tp{Y), respectively, 
where 0o — {(/>) and ipo = (ip) are the spatially averaged 
monolayer compositions, and Scj) and Sip describe the de- 
viations from their average values. In thermal equilib- 
rium, 00 and ijjQ satisfy 



T^<Po + 4'l- pi-4>- AV'o = 0, 



(12) 



Expanding the free energy, Eq. and retaining the 

quadratic order terms in (50 and 5^j, we obtain the Gaus- 
sian free-energy 



d^r 



2B(V2 



+2D{V^5i^f - 2C{y5i:f + ^-{6i;f 



(13) 



where the notations tcf, — + 30q and = + S'^q as 
weU as Eq. (fT2|) are used. 

The static partial structure factor is S'00(q) = 
((50(q)(50(— q)), and because of the radial symmetry, 
S4>4>{q) = S^^{q), where q = |q|. It is given by 



S4,4,{q) 



and similarly for and S^^ 

250(9) 



A2 



Sti^^iq) 



450('ZW('Z)-A2' 



S4>^{q) = S^^{q) 



A 



4ff0(g)5^(g) - A2' 



(14) 



(15) 



(16) 



and 



5^(g) = 2B/-2V + f , 
g4q)^2Dq^-2Cq^+'-^. 



(17) 



Since the structure factors diverge at T^, we see that 
the coupling parameter A effectively shifts the critical 
temperature to lower values. 
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B. Decoupled leaflets (A = 0) 



and the second is the correlation length 



When the two leaflets are decoupled, A = 0, it is suffi- 
cient to present results for only one of the two monolay- 
ers, say the <j) one. From Eq. (ITU) , the decoupled struc- 
ture factor is simply given by 



S4,4,{q) 



1 



(18) 



This function has a peak at = -^Z A/2B where has a 
minimum, indicating that monolayer fluctuations have a 
characteristic length scale describing their spatial modu- 
lations. The peak height is hence given by 



1 



(19) 



and diverges at = /B. In addition, the structure 
factor in Eq. (I18|) decays as ^ for large wavenumbers. 
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FIG. 3. The rescaled (dimensionless) real-space correla- 
tion function, G^^, of a single monolayer in the structured- 
disordered phase as a function of r. The parameters are 
B/e^ = 1 and 7^ = —0.9, and the two resulting characteristic 
lengths are X^/2n ~ 1.45 and ~ 6.32. 

The correlation function G^^{r) — (50(r)50(O)) = 
G^^{r), where r — |r|, is obtained by the inverse 2D 
Fourier transform of Eq. ^TE\\ (see Appendix A for the 
derivation) 



G00(r) 



where _ffg^' (z) is the zeroth order Hankel function of the 
first kind, and "Re" stands for the real part. This corre- 
lation function contains two length scales; the first is the 
modulation periodicity 



32ttB' 



Re 



(1) 



27rr . r 



(20) 



2^ 



1/4 



\/ 1 - 70 ' 



(21) 



^0 



B 



1/4 



(22) 



where 7^ = -A/ ^/e^B. 

In Fig. |3l we plot Eq. (|20)) for certain parameter val- 
ues. The correlation function has exponentially decaying 
oscillations as a function of the distance r, similar to 3D 
microemulsions ^2§\. When — 1 < 7^ < 1, both A^ and 
^0 are finite, and the corresponding phase is called the 
structured-disordered phase. The peak of the structure 
factor occurs at g* > for — 1 < 7^ < 0, whereas 
g* = for < 70 < 1. The line at 7^ = is the Lifshitz 
line. The modulation periodicity A^ diverges for 7^ — >■ 1, 
and the monolayer transforms into a disordered phase for 
70 > 1 . The line of 7^ = 1 is called the disorder line and 
is not a phase transition line. On the other hand, the 
correlation length diverges for 7^ — >■ — 1 (hence, this 
is the critical point), and the ordered phase appears for 
70 < -1 [28^. 

As shown in Fig. Ufa), the phase diagram of a decou- 
pled bilayer is easily obtained by combining the phase 
sequences in terms of the two independent parameters 



70 = —A/^e^B and 7^ = —C/ ^Je.^D. The various 
phases are expressed by the binary combination of or- 
dered (O), structured-disordered (S) and disordered (D) 
phases for each of the leafiets. The two disorder lines are 
shown in the figure by dashed lines, while the Lifshitz 
lines occur at 70 = and 7^ = (not shown). 



C. Coupled leaflets (A / 0) 

When the two monolayers are coupled (A ^ 0), the 
corresponding phase diagram can be obtained by analyz- 
ing the poles in the complex plane of Eqs. (ITH - p^ . One 
example is shown in Fig. SJ^b). The boundary between 
the 00 and SS phases (solid line) is determined by the 
condition that all the structure factors diverge. In the 
DD phase, all the poles are pure imaginary, whereas in 
the SS phase the poles are complex, and at least one real 
pole exists in the 00 phase. 

In Fig. mjb) we see that the asymmetric phases such as 
DO and SD are suppressed as compared with Fig. HJ^a). 
This is because one of the leaflets induces modulations 
in the other leaflet due to the inter-leaflet coupling, A. 
In addition, the 00 and SS phases are expanded as 
compared to Fig. SJ^a) . As mentioned before, increasing 
the coupling A effectively lowers the bilayer temperature. 
The Lifshitz lines (dotted lines) are obtained by plotting 
the main peak positions of Eqs. and ((T5|) . Due to 
the coupling effect, they are tilted as compared to the 
A = case where the lines occur at 70 = and 7^ = 0. 

In Fig. [SJa) and (b), we plot the structure factors of 
the decoupled and coupled bilayers, respectively. As an 
illustration of the coupling effect with A = 0.3, we con- 
sider two different characteristic wavelengths ^ g^. 
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FIG. 4. (a) Phase diagram of a decoupled modulated 
bilayer (A = 0) plotted in the plane of 7^ — —Aj^e^B 
and 7^ — —C / t^D. The phases are labeled using a bi- 
nary combination of two letters representing phases in the 
cj> and monolayers. "O", "S" and "D" stand for the or- 
dered, structured-disordered and disordered phases, respec- 
tively. Solid and dashed lines are the phase transition and 
disorder lines, respectively, (b) Phase diagram of a coupled 
modulated bilayer (A = 1) for = — 2, B = D = 0.5. In 
addition, the two dotted lines are the Lifshitz lines. 



FIG. 5. Bilayer structure factors S^,,/,, S^^, and as a 
function of the wavenumber q for = = 1.5, B = A — 1, 
D = 0.0123, C = 0.1111 = 1/^2,(7; = 3/V2). (a) the 
decoupled case, A = 0; (b) the coupled case with A = 0.3. 



The heights of the two peaks are set to be equal by re- 
quiring £0 = and A^/B = C^/D = 1 (see Eq. (iM ). 
The peak height of S^^ at = 1/ -\/2 is increased (as in 
Fig.EKb)) due to the coupling effect, whereas that of S^^ 
s-t —3/ \pl is almost unchanged compared with the de- 
coupled case. We also plot S^^ given by Eq. ([T6| which 
represents the cross correlation of fluctuations between 
the two monolayers. This quantity is proportional to the 
coupling constant A and its peak position is essentially 
determined by that of S'^^ at = 1/ \/2. 

The peak position of S^^ denoted as for the cou- 
pled case (A = 0.3) is obtained numerically and plotted 
in Fig. inja) as a function of q*^ and q*^ , which are, respec- 



tively, the peak positions of and for the decou- 
pled case. We find that for the coupled case, the value 
of gj,^ is almost equal to that of the smaller of (gj, g^). 
In Fig. [6l^b) , the peak height of S^^ at (f^.^ is plotted as 
a function of the coupling parameter A and the ratio of 
the two characteristic wavenumbers g^/gj for the spe- 
cific case of gj = 1/a/2, as the parameters are fixed to 
B — A — \ and jD = 1. The peak height increases 
as A is increased because the temperature of the coupled 
bilayer is effectively lowered. To see this clearly, we have 
plotted in Fig.lHJb) the phase-transition line (dashed line) 
separating the structured-disordered and ordered phases. 
The peak height of S,^^ also increases as q^/q*^ = V^^i/, 
approaches unity, which is the case where S'^^ and S^^ 
completely overlap. 
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FIG. 6. (a) Contour plot of the peak position 



of 



where and are the kinetic coefScients taken to be 
constants, and t'i(r, represent Gaussian white noise, 
satisfying {yi[v.t)) — and 

{y,(v,t)v,(v\t')) = -%L,v25(r - r')<5(i - t'), (24) 

where z, j — (j),ip and (• • • ) now indicates the average over 
space and time. 

The Fourier transform of fluctuations in both space 
and time is 

6<j){<i,uj)= J d^rdtS(t>{r,t)e-'^'^"-'^'\ (25) 

and a similar Fourier transform is used for dip and 
j/j. The dynamic structure factors such as S^^{q,uj) — 
{S(f>{q,uj)S(f){~q,—uj)) are given by 



Si] (q, w) = 



aij{q) 



;2 + [c._(q)]2' 



(26) 



where the explicit expressions of aij{q) and Pij{q) are 
given in Table |T] and the characteristic frequencies uj± (q) 
are: 



[^±{q)? = l 



'^l{q) + ^l{q) 



Luliq) = ALyMq)]^ + L^L^q^k\ 
uliq) = ALlq%^{q)f + L^L^q^K\ 
^4>i>iq) = '2-q^[L,t>9<t>{q) + L^g^{q)], 



(27) 

(28) 
(29) 
(30) 



the str ucture factor S^^, as function of = ^ AjlB and and qa. and q,;, are defined in Eg. (flT)) . In the more special 

= /2D. The parameters are A = 0.3, — e^p = 1.5, case of = - 
and A^/B — C^/D. The dotted line represents — gj. 

(b) Contour plot of the peak height of the structure factor ^ 
S^Tp{q*^^) as a function of A and %/2qJ. For the choice of '-^±{q) — Lq 

parameters B = A = 1 and /D = 1, ~ l/\/2, the 

two characteristic wavenumbers satisfy gj/g^ = \/2gJ. The ip ^/ {g^{q) - g^{q)Y + . (31) 

dashed line is the phase-transition line. 



L, w± reduce to a simpler form 

94>{q) + 9^{q) 
T i94>{q) - 94:{qW + 



D. Dynamic structure factors TABLE I. Expressions for aij(g) and Pij{q) 



ajjil) Pij (g) 



The dynamical fluctuations in composition, S(j)(r, t) 
and t), depend on time t and are now considered for 2L^q'^{uj^ — uj\) 2L^q'^{w^ — w^) 
coupled modulated bilayers. Since both S(j) and S-tp are '^'i"!' ^ u?~^u? '^'^'^ ^ a?~^a? 
conserved quantities for each monolayer, the time evolu- 
tion equations are given by 07-2/2 2 ^ r,r2/2 2x 

a^V — 2 2 Pi'i' — 2 1 

U)_ — OJj^ ui_ — U)^ 

d5(l>{Y, t) _ ^2 &Fg 

^|^=i.V^^ + .,(r,0, (23) 
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FIG. 7. Decay time of concentration fluctuations as a func- 
tion of the wavenumber q for — = 1.5, B — A = 1, 
D = 0.0123, C = 0.1111 (resulting in q*^ = 1/^2 and 
q^ = 3/-\/2), and L = 1. (a) The decoupled case, A = 0; 
(b) the coupled case with A = 0.3. 



The mtermediate structure factors Sij{q, t) depend ex- 
plicitly on time and are obtained by taking the inverse 
Fourier transform in t: 



and from Eq. (\26^ . we obtain 



2io_iq)' 



(32) 



(33) 



Hence the decay of concentration fluctuations is de- 
scribed by a sum of two exponentials with two character- 
istic times 1/Lj±{q). However, when the two monolayers 
are decoupled (A = 0), the S^^ and S^^ structure factors 
decay with a single exponential characterized by a decay 
rate and (jJ^ (see Eqs. (1^ and dH])), respectively. 

For the decoupled case (A = 0), the decay times l/w^ 
and I/uItp are plotted as a function of q in Fig. [71[a), 



with the same parameters as those in Fig. [5] The plots 
show a shoulder reflecting the characteristic structure at 
wavenumbers = l/\/2 and = 3/-\/2. Notice that 
the larger the initial length scale, the longer the decay 
time. For the coupled case with A = 0.3, we plot l/aj± 
in Fig. [7fb). Due to the coupling effect, the two decay 
times split into a larger and smaller one, > 
The larger one, exhibits two shoulders, while the 

smaller one, has a shoulder between the two char- 

acteristic wavenumbers. The coupling affects the de- 
cay time of the structure corresponding to the smaller 
wavenumber (larger length), similar to the effect seen in 
Fig. E] for the static structure factor. 



IV. COUPLED MODULATED PHASES WITH 
DIFFERENT PERIODICITIES q^ / q^ 

In our previous paper j29| . we investigated the phase 
behavior of coupled modulated bilayers by using the free 
energy, Eq. (|lip . below T^. When the two monolayers 
have the same preferred periodicity, q*^ = q^, we ob- 
tained the mean-field phase diagram exhibiting various 
combinations of modulated structures such as stripe (S) 
and hexagonal (H) phases. In some cases, the periodic 
structure in one of the monolayers induces a similar mod- 
ulation in the second monolayer. Moreover, the region of 
the induced modulated phase expands as the coupling 
parameter A becomes larger. 

However, when the preferred periodicities in the two 
leaflets are different, q'^ Q^, it is difficult to obtain the 
phase diagram because the free energy densities cannot 
be obtained analytically. We then have to rely on numer- 
ical simulations [2^ Hlj to solve the time evolution of the 
two coupled order-parameters as explained below. 

For the dynamics of coupled modulated bilayers below 
Tc, we use the following time evolution equations [2^ 



dt 
dt 



(34) 



where the bilayer free-energy, i^b, is given by Eq. 
Below Tc the effect of thermal fluctuations is less impor- 
tant and the noise terms have been omitted in the above 
equations. Hereafter, the kinetic coefficients and 
are set to unity for simplicity. 

We solve numerically the above 2D equations using pe- 
riodic boundary conditions. Each run starts from a ho- 
mogeneous state with a small random noise around the 
average compositions and ipQ. Time is measured in 
discrete time steps, and t = 5000 corresponds to a well 
equilibrated system. In all our simulations, the reduced 
temperature parameters are fixed to be — = 0.8. 
The characteristic wavenumber in the <j) monolayer is 
fixed as q*^ = l/v^ by setting B = A ^ 1, while the 
periodicity in the ip monolayer q^ — \J C /2D is varied 
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FIG. 8. Time evolution of spatially modulated patterns in a coupled bilayer, which consists of two stripe monolayers with 
different periodicities, q'^ 7^ gj. Real-space patterns of (j>, ip, <^ + V') as well as the Fourier transform, ip{c\), are presented. 

The time steps are: (a) t — 25, (b) t = 60, (c) t = 250, (d) t — 5000. The parameters are tpo = ipo = 0, = = 0.8, 
B — A = 1, D = 0.0123, C = 0.1111, A = 0.3, yielding q'^/q*, = 3. In order to emphasize the color contrast, the gray color 
code of the real-space patterns is chosen to vary between —1.7 and 1.7. 



and the condition jD = 1 is used as before. In the 
following, we shall consider only the two coupled stripe 
phases for 0o = V'o = 0- However, the combination of 
stripe and hexagonal phases leads to a rich variety of 
complex patterns [41]. 

When the characteristic wavenumbers of the two de- 
coupled monolayers are different, q*^ ^ gj, the two cou- 
pled modulated structures cannot apparently match each 
other. The frustration between the two different peri- 
odicities is due to the inter-leaflet coupling and affects 
their morphologies. In Fig. [Sj we show one example of 
the time evolution of a coupled micro-phase separation 
for 'Z;/'?; = 3 (g; = Vv^ and q% = ?,/V2) and A = 0.3. 
The spatial patterns of 0, "0, <j) + ipj 4'~4' f^^id the Fourier 
transformed pattern, '(/'(ct)j presented for time steps 
oi t = 25,60,250 and 5000. Starting from the isotropic 
state, the ip monolayer forms first stripes (a). Then, as 
the (j) monolayer starts to segregate at larger t, it simul- 
taneously chops the tp stripes into smaller sections (b), 
and the ip monolayer transforms into a finger-like pat- 
terns (c). Reconnection of the (j) stripes takes place after 
a long time of annealing, and a pattern of alternating fin- 



gers in the tp monolayer is finally obtained (d) . From the 
time evolution of the 2D Fourier patterns of the ip mono- 
layer, it is apparent that the intermediate structures are 
characterized by the two length scales of ratio 1:3. 

In Fig. ini we show the spatially modulated patterns at 
t ~ 5000 of the two coupled (j> and ip monolayers with dif- 
ferent periodicities (q^ = and = 3/\/2 as before) 
for different values of the coupling parameter A = 0.1, 0.3 
and 0.5. Notice that Fig. EJb) is the same as Fig. M^d). 
In the weak-coupling case (A — 0.1), the two monolayers 
exhibit two independent stripe morphologies with charac- 
teristic periodicities (called the "independent" morphol- 
ogy). Here the two stripes essentially do not affect each 
other. As the coupling constant is increased (A = 0.3), 
stripes with a finger-like structure appear in the ip mono- 
layer, while the stripe morphology in the (p monolayer 
is almost unaffected (called the "intermediate" morphol- 
ogy). In the Fourier transformed pattern of Tp, we clearly 
see that the structures with two different characteristic 
wavenumbers are coexisting. This result is also in accord 
with the properties of the static structure factors, Sij{q), 
having two different characteristic lengths, as shown in 




FIG. 9. Equilibrated patterns (time step t — 5000) of two coupled bilayers consisting of stripe monolayers with different 
periodicities. All parameters beside A are the same as in Fig. |8l The A parameter values are: (a) A = 0.1; (b) A = 0.3; (c) 
A = 0.5. In order to emphasize the color contrast, the color code of the real-space patterns is chosen to vary between —1.7 and 
1.7. 



Fig. [5l For a larger coupling parameter (A — 0.5), very 
similar patterns are obtained for (f> and "0, and almost 
coincide with one another (called the "coincident" mor- 
phology). It should be noted that the structure with 
the larger wavelength dominates when the coupling is 
large enough, and the sequence of morphological changes 
shown in Fig. |9l as A increases, is rather typical. 

In order to quantify the three morphologies (indepen- 
dent, intermediate, coincident), we calculate the spatial 
average of the product of the two compositions 

{H) = ^j A\<l,{v,t),p{v,t), (35) 

where A is the total system area. In Fig. [101 we plot 
{(pip) t — 1000 (sufficient for equilibration) for vari- 
ous combinations of A and q^/q"^ — V^q^ (for the case 

— l/V^). The morphology of the obtained patterns 
is marked on the figure by circles, crosses and squares. 
The values of {(pip) are small for "independent" (circles) 
structures, while they become larger for the "intermedi- 
ate" (crosses) and "coincident" (squares) morphologies. 
Although the morphological changes are gradual and do 
not represent a sharp transition, the intermediate pat- 
terns appear roughly for 0.05 < {(pip) < 0.2. The re- 
gion of intermediate structure expands as q^/q*p ~ V^q^ 
increases, and the patterns coincide for A = 0.5. We 
note that although the morphology cannot be solely de- 
termined by the quantity {(pip)^ the behavior of {(pip) is 



similar to that of the peak values of the cross correlation, 
Sui'tu'^-' presented in Fig.^b). 

For the quantitative argument concerning the micro- 
phase separation dynamics in each monolayer, we have 
also calculated the two self quantities: 

{(P') = ^j d\(t>\v,t), (36) 

and 

{V) = ^j A\ip\v,t). (37) 

In Fig. [11] we plot the square root of these quantities as 
a function of with the same parameters as those used 
in Fig. [9l In all studied cases, the modulation of the ip 
monolayer having a larger wavenumber (g^ > gp, grows 
faster than that of the (p monolayer. We also remark 
that the structure formation in the (p monolayer is accel- 
erated for larger A, whereas that of the ip monolayer is 
almost unchanged. According to the linear stability anal- 
ysis of Eq. (|34|) . the initial growth rates of the unstable 
modes are essentially equivalent to the decay rates uj± of 
the concentration fluctuations given in Eq. (j3ip . This is 
consistent with Fig. [Tj where the characteristic growth 
time (l/cj) for larger q is smaller than that for smaller 
q. The growth rates increases with A (see Fig. [TlTa)) 
because the coupling effectively reduces the temperature 
and enhances the phase transition. The faster the decay 
of concentration fluctuations, the faster the structure for- 
mation. 
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FIG. 10. Contour plot of ((^^/') at t = 1000 as a function 
of A and \/2gJ, for the specific case of = l/\/2. (a), (b), 
(c) correspond to A = 0.1,0.3 and 0.5, respectively, as given 
in Fig. [9] Circles, crosses, squares correspond, respectively, 
to independent, intermediate, and coincident structures (see 
text). 



V. DISCUSSION AND FINAL REMARKS 

In this paper we present a model for coupled modu- 
lated lipid bilayers. We start by considering a monolayer 
consisting of a mixture of saturated and hybrid lipids, 
and propose a phenomenological model that includes a 
coupling between the lipid composition and a 2D vecto- 
rial field. This coupling arises from the line-active nature 
of the hybrid lipid, which adjusts its tail orientation in or- 
der to reduce the line tension. Minimization of the mono- 
layer free-energy with respect to the vectorial field yields 
a 2D microemulsion exhibiting modulated phases. The 
characteristic wavelength of modulation is determined by 
the monolayer coefficients A and B, Eq. ([8]), reflecting 
molecular properties of lipid mixture. We then construct 
a model for lipid bilayers comprised of two modulated 
monolayers that influence each other through an inter- 
leaflet coupling. 

Based on the model, we study concentration fluctua- 
tions of bilayers above Tc, and calculate their static struc- 
ture factors. The calculated phase diagram for coupled 
bilayers shows that the extent of ordered and structured- 
disordered phases become larger as compared to the de- 
coupled case. When the two monolayers with differ- 
ent preferred wavenumbers ^ are coupled, (say 
14> ^ 90 ); the peak height of S^^ occurring at smaller 
q-numbers becomes larger as compared to the decoupled 
case, whereas the peak height of S,p,p occurring at larger 
q-numbers almost does not change. Namely, the inter- 
leaflet coupling strongly affects the compositional modu- 
lation in each monolayer. Furthermore, the inter-lcaflct 
coupling has a clear signature on the cross structure fac- 
tor, S^^, as well as on the dynamics of concentration 




FIG. 11. Time evolution of the modulation amplitudes of the 
two monolayers: (a) ^ {(p^) and (b) ^ (tp^) . The parameters 
are = = 0.8, <f>o = i>o = 0, B = A = 1, D = 0.0123, 
C — 0.1111 yielding q'^/q^ — 3. The different lines represent 
A = 0.1 (solid), A = 0.3 (dashed), and A = 0.5 (dotted). 



fluctuations. By calculating the intermediate structure 
factor, S{q,t), we show that concentration fluctuations 
generally exhibit a double exponential decay with two 
decay rates, uj±. One of the decay times (1/w-i-) exhibits 
two shoulders at wavenumbers describing the monolayer 
compositional modulations. 

For membranes below Tc, we studied the micro-phase 
separation of a coupled modulated bilayer. When the 
two monolayers have different modulations, 7^ g^, we 
obtained numerically a variety of complex patterns. The 
initial growth rates of the unstable modes are identical 
to the decay rates of the concentration fluctuations. 

As mentioned in Sec. I, the special character of concen- 
tration fluctuations in our model may explain the flnite- 
size domains ("rafts") in biological membranes. For or- 
dinary binary mixtures above Tc, the only length scale 
of the disorder phase is determined by the correlation 
length, and close to Tc, this length becomes large. But 
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within our 2D microemulsion model, there is another 
length scale characterizing the modulations as given by 
Eq. (|21|) . This second length scale may also explain the 
finite-size domains observed in some experiments as a 
result of micro-phase separation in the low temperature 
phase. 

A related model based on a microemulsion picture was 
recently proposed by Schick [4^, who considered a cou- 
pling between curvature and compositional asymmetry 
between the two leaflets, resulting in a 2D microemul- 
sion. Although Schick's model as well as ours share the 
microemulsion viewpoint, the origin of the physical mech- 
anism is different. In Schick's model the coupling be- 
tween composition asymmetry and curvature gives rise 
to domains with different spontaneous curvature. Our 
model, on the other hand, assumes that for flat mono- 
layers, the line-active nature of hybrid lipids is solely re- 
sponsible for the microemulsion formation (see Eq. ([3l)). 
Hence, our proposed physical mechanism for the bilayer 
coupling, as discussed in Sec. IIIBl is different. One of 
the consequences of our A > coupling is that the do- 
mains residing on the two leaflets are in register, unlike 
the prediction of Ref. [i^ . 

The present work is concerned with the analysis of 
equilibrium properties and relaxation dynamics towards 
the equilibrium state. The existence of finite-size do- 
mains may also be explained by non-equilibrium lipid 
transport between the cell interior and the membrane. 
Such a mechanism was considered j43l - l46l | through a cou- 
pling between the membrane and an outer reservoir of 
lipid or cholesterol. Similarly to our model, these works 
explain the appearance of the finite-size domains as a 
result of micro-phase separation. 

Asymmetry between lipid composition in the inner and 
outer leaflets of biomembranes has a very deep signifi- 
cance and is closely related to the cell biological func- 
tions. For instance, the breakdown of such compositional 
asymmetry is related to programmed cell death (apopto- 
sis) [l^]. For living cells, this asymmetry is maintained by 
an enzyme called "fiippase" that actively transports the 
lipids between the two leaflets [11] . The half-time of lipid 
composition due to flip-flop motion is measured using 
time resolved small angle neutron scattering (SANS) j49| . 
and is estimated to be several hours at physiological con- 
ditions. 

Although membranal signal transduction is important 
for various biological functions, its dynamical properties 
are still not so well understood. As our theory offers an 
explanation for the size and dynamics of lipid domains, 
we hope that it and similar models will contribute in the 
future towards the understanding of functional processes 
in biomembranes. 
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Appendix A: Derivation of G^^(r) 

In this Appendix we present the derivation of G^^ of 
Eq. (PO)) . The real-space correlation function is given by 
the 2D inverse Fourier transform of Eq. (fT8)) : 



G^^{r) 



d^g 
(2^)2 



AttB Jo ^ q^-{A/B)q^ + e^/AB' 



(Al) 



where Jo{qr) is the zeroth-order Bessel function. We 

use the relation Jo{qr) — [HQ^\qr) + H'^\qr)\/2 where 

H^^\qr) and H^^\qr) are the zeroth-order Hankel func- 
tions of the first and second kind, respectively. Then, the 
integral in Eq. (|A1|) is written as / = (/i -|- J2)/2 where 



with i = 1, 2. 



dq- 



qH^'>iqr) 



{A/B)q^ + €^/AB' 



(A2) 




FIG. 12. Four poles in the complex plane at z = Zj {j — 
1,2,3 and 4). The two close integral paths are indicated by 
the arrows in the 1st and 4th quadrants. 

In order to evaluate the above integral, we performed 
the integration in the complex plane by replacing q with 
the complex variable z = x -\- iy. The integrand has four 
poles at 



(e^/i?)^/^ 



(±v/l-70±«v/l + 70) 



±- 



27r 



±i- 



(A3) 
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with j = 1,2,3,4 and and have been defined in 
Eqs. (PT|) and (P^ . These poles are located in the quad- 
rants 1, 2, 3 and 4, off the real and imaginary axes. For 
the integral /i, we integrate along a quarter-circle con- 
tour in the first quadrant in an anti-clockwise direction. 
The contour radius is taken to infinity (Fig. I12p. and us- 
ing the residue theorem we obtain 



oo 



da; ■ 



^ 



+ / dy 



■4 - {A/B)x^ + e^/AB 
.JlS^rJ ,,;, =2^^Res|.., (A4) 



where "Res" stands for the residue. Similarly, for the in- 
tegral /2, we integrate along the contour of the quarter- 
circle of infinite radius in the fourth quadrant in the 



clockwise direction: 

xHI^^ (xr) 



dx 



x'^ - {A/B)x^ + e^/AB 



dy 



iyH^^\iyr) 



y4 + {A/B)y^ + e^/4B 



-27riResU^. (A5) 



Combining Eqs. (jA4p and (|A5|) . and further using the 
relation Hq^\—z) — —H^\z), we obtain 



1 = 



H^^\z^r)+H^'>{zir) . (A6) 



Finally, using and from Eq. (jA3p . we obtain 



" A^ 



(A7) 



In the above expressions, we have used the relation 
of z. 



h'^\z) = hI^\z), where z is the complex conjugate 
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